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Abstract 

To  achieve  a  better  understanding  of  the  degradation  phenomena  of  polymer  electrolyte  membrane  fuel  cells  (PEMFCs),  it  is  imperative  to 
understand  the  mechanism  of  microstructure  changes  in  the  catalyst  layer.  To  this  end,  a  rate-dependent  isotropic  plasticity  model  with  temperature- 
and  humidity-dependent  material  properties  is  proposed  to  describe  the  viscoplasticity  of  the  catalyst  layer  components.  To  understand  the 
mechanism  of  such  changes  caused  by  the  cycling  of  start-up  and  shutdown  during  the  operation  of  a  PEMFC,  the  material  model,  combined  with 
the  cohesive  zone  model  and  the  contact  model,  is  solved  using  the  finite  element  method.  The  cohesive  zone  model  and  the  frictional  contact 
model  are  used  to  describe  the  evolution  of  interfaces  between  the  protonic  and  the  electronic  conducting  phases.  Numerical  simulation,  based 
on  the  representation  of  the  micro  structure  in  the  catalyst  layers,  shows  that  there  is  competition  between  crack  initiations  in  the  bulk  material  of 
the  protonic  phase  and  delamination  between  different  phases.  This  competition  plays  an  important  role  in  micro  structure  changes  in  the  catalyst 
layers.  The  reduction  of  connectivity  between  the  protonic  and  the  electronic  conducting  phases,  which  may  explain  the  decrease  of  performance 
after  certain  duty  cycles,  could  be  contributed  to  by  cracking  or  delamination. 
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1.  Introduction 

The  durability  of  polymer  electrolyte  membrane  fuel  cells 
(PEMFCs)  is  one  of  the  most  important  issues  in  the  push 
to  successfully  commercialize  broad  stationary  and  transporta¬ 
tion  energy  applications.  In  fact,  the  lifetime  of  PEMFCs  under 
steady  operating  conditions  is  about  three  to  five  times  longer 
than  under  dynamic  operating  conditions  [1,2].  For  example, 
the  lifetime  of  PEMFC  stacks  has  reached  over  10,000  h  in  sta¬ 
tionary  applications;  by  contrast,  state-of-the-art  PEMFC  stacks 
can  be  run  for  only  2000  h  in  automotive  applications  [3-7]. 
Therefore,  it  is  urgent  to  better  understand  why  the  lifetime  in 
stationary  and  dynamic  operating  environments  is  so  different 
and  how  to  overcome  it. 
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Recently,  experimental  studies  have  been  undertaken  to 
understand  the  mechanism  of  degradation  and  to  characterize  the 
aging  phenomena  as  functions  of  time.  Many  researchers  [3-6] 
have  indeed  observed  micro  structure  changes  when  comparing 
aging  membrane  electrode  assemblies  (MEAs)  with  fresh  ones. 
These  changes  or  failures  in  MEAs  result  in  a  gradual  reduc¬ 
tion  of  ionic  conductivity,  an  increase  in  total  cell  resistance, 
a  reduction  of  active  sites,  a  reduction  of  voltage,  and  a  loss 
of  output  power  [4].  As  crucial  parts  of  MEAs,  catalyst  layers 
(CLs)  are  examined  post-mortem  through  transmission  electron 
microscopy  (TEM)  and  scanning  electron  microscopy  (SEM). 
Many  researchers  believe  that  the  microstructure  changes  in 
CLs  can  be  seen  as  factors  in  fuel  cell  performance  reduction 
[8,9].  According  to  experiment  results,  CLs  show  the  follow¬ 
ing  microstructure  changes:  cracking  or  delamination,  loss  of 
carbon- supported  catalyst  clusters,  dissolution  of  the  electrolyte 
(Nafion  ionomer),  catalyst  particle  migration,  catalyst  ripening, 
carbon  coarsening,  and  so  on  [5,10].  These  phenomena  are  rela¬ 
tively  obvious,  especially  under  dynamic  operating  conditions. 
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So  the  microstructure,  especially  in  CLs,  is  arguably  a  critical 
point  for  improving  the  durability  of  PEMFCs  and  more  and 
more  researchers  are  focusing  on  elucidating  the  mechanism  of 
microstructure  changes  in  CLs  and  their  effects  on  performance. 

Microstructure  changes  in  CLs  may  occur  in  several  ways, 
such  as  chemical  degradation  of  the  ionic  conducting  parts 
or  mechanical  failure  [3-5].  Although  chemical  degradation 
is  a  key  factor  in  microstructure  changes,  it  has  been  sug¬ 
gested  that  mechanical  damage  also  plays  an  important  role 
[11,12].  Lor  example,  debonding  between  the  electrolyte  and  the 
carbon-catalyst  agglomerate  directly  results  in  a  loss  of  carbon- 
supported  catalyst  clusters,  while  the  breaking  of  the  electrolyte 
network  contributes  to  the  dissolution  of  the  proton-conducting 
network  [5].  Karlsson  et  al.  [12,13]  used  a  model  to  simulate  the 
mechanical  analysis  of  the  MEA  scale.  As  a  result  of  the  start¬ 
up  and  shutdown  of  the  fuel  cells,  humidity  and  temperature 
fluctuate  and  influence  the  mechanical  behaviour  of  different 
MEA  components.  Karlsson  et  al.  found  that  critical  residual 
stresses  accumulate  and  lead  to  mechanical  fatigue  in  the  mem¬ 
branes  after  a  hydrothermal  cycle.  However,  the  mechanism 
of  microstructure  changes  in  CLs  is  still  unclear.  Mechani¬ 
cal  damage  in  the  CL  can  appear  as  flaws  or  mud-cracks  [4], 
delamination  between  the  carbon-catalyst  agglomerate  and  the 
electrolyte  [10],  or  corrosion  of  carbon  [6].  Cycled  dynamic 
operating  conditions  should  result  in  such  mechanical  dam¬ 
age  or  degradation,  and  this  will  greatly  affect  the  decline  in 
performance. 

As  discussed  above,  although  many  experiments  show  that 
microstructure  changes  occur  in  CLs  after  long,  cycled  operation 
of  PEMFCs,  so  far  there  is  no  theoretical  or  numerical  model 
to  describe  either  the  chemical  or  the  mechanical  aspects  of  this 
evolution  in  structure.  The  purpose  of  this  work  is,  first,  to  estab¬ 
lish  a  theoretical  framework  and  a  numerical  model  to  investigate 
microstructure  changes  in  CLs,  and,  second,  to  reveal  the  chang¬ 
ing  mechanism  and  to  understand  its  effects  on  performance.  As 
a  first  step,  in  this  paper  we  focus  only  on  the  mechanical  aspects 
of  microstructure  changes  in  CLs.  The  effects  of  such  changes 
on  CL  performance  and  the  modeling  of  the  changes  caused 
by  chemical  property  variation  will  be  discussed  in  subsequent 
papers. 

In  this  paper,  the  primary  methodology  framework  proposed 
is  the  finite  element  method  (FEM);  the  cohesive  zone  model 
and  the  contact  model  are  used  to  describe  responses  of  bulk 
material,  interfacial  behaviours  between  Nation  and  the  carbon- 
catalyst  agglomerate  and  contact  phenomena  between  solids  in 
CLs,  respectively.  The  temperature-,  humidity-  and  strain  rate- 
dependent  material  models  for  Nation  are  based  on  a  literature 
review  and  empirical  results.  A  simulation,  based  on  a  represen¬ 
tation  of  CLs,  is  performed  and  shows  that  the  failure  of  Nation 
or  the  delamination  between  Nation  and  the  C/Pt  agglomerate 
can  be  caused  by  hydrothermal  cycles  even  in  the  absence  of  a 
chemical  reaction.  The  mesh  dependency  of  the  present  model 
is  examined  through  the  simulation  of  different  mesh  densities. 
Moreover,  a  parametric  study  is  done  to  investigate  the  effects  of 
different  frequencies  of  start-up  and  shutdown  of  the  PEMFC.  A 
possible  indication  of  performance  decline  is  observed,  based  on 
the  phase  connection  between  Nation  and  the  C/Pt  agglomerate. 


Finally,  through  a  mechanical  analysis,  a  basic  understanding  of 
the  mechanism  of  microstructure  changes  in  CLs  is  achieved. 

2.  Basic  framework 

The  main  focus  of  this  paper  is  the  mechanical  mechanism  of 
the  evolution  of  CL  microstructure.  When  a  PEMFC  is  starting 
up  or  shutting  down,  the  humidity  and  temperature  in  its  CLs 
increases  or  decreases,  respectively.  When  the  cell  is  operating 
steadily,  the  humidity  and  temperature  are  relatively  stable.  The 
change  in  relative  humidity  (RH)  or  temperature  ( T )  is  illustrated 
in  Fig.  1(a).  In  order  to  simulate  microstructure  changes  caused 
by  the  duty  cycles  of  fuel  cells,  we  ignore  the  creep/relaxation 
processes  occurring  in  part  B  and  take  into  account  only  the 
relatively  fast  material  fatigue  process  occurring  in  parts  A  and 
C.  So  only  the  start-up  and  shutdown  processes  of  PEMFCs 
are  considered,  and  the  mechanical  analysis  is  based  on  these 
cycles  of  RH  and  T ,  shown  schematically  in  Fig.  1(b)  and  (c). 
The  periods  of  RH  and  T  cycles  are  £rh  and  tj,  respectively. 
Based  on  this  assumption,  one  mechanical  model  is  proposed  to 
investigate  the  effect  of  RH  and  T  cycles,  simulating  the  start-up 
and  shutdown  of  the  fuel  cell. 

2.7.  Overview 

There  are  three  different  phases  in  the  CLs  of  a  PEMFC: 
electrolyte  (Nation),  pore,  and  carbon-catalyst  (Pt,  etc.)  agglom¬ 
erate.  As  indicated  in  Fig.  2,  TEM  is  utilized  to  examine  the 
microstructure  of  CLs  [5]. 

The  pore  phase  is  considered  void  in  the  solid.  The  solid 
phases  (Nation  and  C/Pt  agglomerate)  are  meshed  in  the  domains 
for  the  sake  of  the  FEM.  The  interfaces  between  any  two  of 
the  three  phases  gradually  change  during  the  operation  of  the 
PEMFC  when  a  cycled  load  is  applied.  Generally,  debonding  or 
delamination  can  occur  between  the  Nation  and  the  C/Pt  agglom¬ 
erate  phases,  so  the  cohesive  zone  model  (CZM)  is  used  to  mimic 
the  evolution  of  the  interfaces  between  these  two  solid  phases, 
indicated  by  the  dashed  line  in  Fig.  2.  The  boundaries  between 
pore  and  solid  change  because  of  the  deformation  of  the  solid. 
So  the  contact  model  is  applied  on  these  interfaces  (indicated  by 
the  solid  lines  in  Fig.  2). 

2.2.  Mechanical  analysis  model  for  bulk  material 

We  start  with  the  FEM;  detailed  procedures  can  be  found 
in  Ref.  [14].  The  material  experiments  already  show  that  the 
mechanical  properties  of  CLs  exhibit  a  plastic  or  inelastic 
behaviour  during  PEMFC  operation.  This  plastic  behaviour  of 
the  materials  has  to  be  considered,  as  does  their  elasticity,  in 
order  to  simulate  the  changes  in  microstructure.  For  this  pur¬ 
pose,  the  inelastic  response,  incompressible  plastic  deformation, 
is  assumed  in  our  model  and  plasticity  theories  are  used  to  char¬ 
acterize  the  elastoplastic  response  of  the  CL  materials.  The  three 
key  factors  in  the  plasticity  theory  of  CL  materials  are  the  yield 
criterion,  the  flow  rule,  and  the  hardening  rule. 

First,  the  von  Mises  yield  criterion  is  adopted  as  the  yield 
criterion  of  the  CL  materials.  It  states  that  yield  occurs  when  the 
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Fig.  1.  (a)  Schematic  cycled  changes  of  humidity  (RH)  and  temperature  ( T)  in  CLs.  Parts  A  and  C  represent  the  start-up  and  shutdown  processes  of  fuel  cells, 
respectively;  part  B  represents  the  steady  operation  of  fuel  cells;  (b)  the  simplified  RH  duty  cycles,  in  which  the  period  is  £rh;  (c)  the  simplified  T  duty  cycles,  in 
which  the  period  is  tj,  and  the  phase  difference  between  RH  and  T  cycles  is  to . 


von  Mises  yield  function  (J2-flow  theory  [15])  approaches  zero. 
The  von  Mises  yield  function  is  given  as 


m  = 


-SS  -  ffp, 


(1) 


where  cr  is  the  stress  tensor,  crp  the  yield  stress  of  the  CL 
materials,  and  S  is  the  deviatoric  stress.  Therefore,  according 
to  the  von  Mises  yield  criterion,  yield  occurs  when  /(cr)  =  0. 
For/(tr)<0,  the  CL  material  exhibits  an  elastic  behaviour  and 
deforms  elastically. 


Fig.  2.  TEM  image  of  different  phases  in  a  CL,  where  the  dashed  line  is  the 
boundary  between  Nahon  and  the  C/Pt  agglomerate,  and  the  solid  lines  are  the 
boundaries  between  the  pore  and  the  other  two  components. 


Second,  the  isotropic  von  Mises  flow  rule  and  the  isotropic 
(work)  hardening  rule  are  used  here.  Flow  theory  predicts  the 
direction  of  plastic  straining  increment  tensor  and  the  hardening 
rule  describes  the  changing  of  the  yield  surface  with  progressive 
yielding,  so  that  the  conditions  (i.e.,  stress  states)  for  subsequent 
yielding  can  be  established.  Therefore,  the  plastic  strain  epl  must 
be  the  strain  that  occurs  in  a  direction  normal  to  the  yield  surface 
and  proportional  to  the  deviatoric  stress  tensor, 

depl  =  ScIa.  (2) 

where  dX  is  a  scalar  proportionality  factor. 

To  introduce  the  effects  of  temperature  and  humidity  on  the 
behaviour  of  CL  materials,  an  uncoupled  theory  is  assumed 
here  in  which  the  additional  temperature  changes  brought  by 
the  strain  are  ignored.  So  the  total  strain  tensor  can  be  given  as 
follows  [13]: 

e  =  ee  +  epl  +  eS  +  eT ,  (3) 

where  ee  is  the  elastic  strain,  epl  is  the  plastic  strain,  and  8s and 
eT  are  the  strains  induced  by  swelling  and  temperature,  respec¬ 
tively.  Let  To  be  a  reference  temperature  and  T  the  current 
temperature;  then  the  thermal  strains  in  Eq.  (3)  resulting  from  a 
change  in  temperature  of  an  unconstrained  isotropic  volume  are 
given  by  [13] 

eT  =  aI(T  -  T0), 

where  a  is  the  coefficient  of  thermal  expansion  (CTE)  and  I 
is  the  identity  tensor.  Similarly,  the  swelling  strains  caused  by 
moisture  uptake  are  given  by  [13] 

es  =  pl( RH  -  RH0), 

where  RH  is  the  relative  humidity  and  RHo  is  the  reference  value 
for  RH.  is  the  coefficient  of  moisture  dependent  expansion 
(CME)  and  I  is  the  identity  tensor. 
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So  far,  an  isotropic  plasticity  model  with  temperature-  and 
humidity-dependent  material  properties  is  proposed  for  the  CL 
materials.  However,  under  practical  operating  conditions,  dif¬ 
ferent  frequencies  of  the  start-up  and  shutdown  of  fuel  cells  will 
have  dramatically  different  effects  on  the  mechanical  response 
of  CL  materials.  The  isotropic  plasticity  model  discussed  above 
does  not  include  the  effect  of  the  strain  rate,  which  is  an  indi¬ 
cation  of  this  frequency.  Many  experiments  on  the  mechanical 
responses  of  polymer  have  been  done  to  clarify  its  viscoplas¬ 
ticity.  Thus  the  isotropic  plasticity  model  should  be  modified  to 
a  rate-dependent  plasticity  or  viscoplasticity  one.  The  deforma¬ 
tion  of  materials  is  assumed  to  develop  as  a  function  of  the  strain 
rate  (or  time).  Many  models  are  available  to  describe  viscoplas¬ 
ticity,  such  as  the  Peirce  model  [16]  and  the  Perzyna  model  [17]. 
In  this  paper  we  adopt  Peirce’s  rate-dependent  plasticity  model, 


£Pi 

1  +  — 
Y 


tfpO, 


(4) 


where  ap  is  the  material  yield  stress,  £pl  the  equivalent  plastic 
strain  rate,  m  the  strain  rate  hardening  parameter,  y  the  material 
viscosity  parameter,  and  apo  is  the  static  yield  stress  of  the  mate¬ 
rial,  which  means  the  yield  stress  under  the  quasi-static  loading 
condition.  Here,  two  parameters,  m  and  y,  must  be  estimated 
from  experimental  data;  this  will  be  discussed  in  Section  3.  The 
reason  to  choose  this  model  over  the  others  is  the  fact  that,  for 
small  values  of  m,  this  model  shows  much  better  convergence. 
The  von  Mises  yield  function  (Eq.  (1))  is  also  rate-dependent 
because  of  viscoplasticity  (Eq.  (4));  that  is, 


3 

gP1 

\  — ss  — 

i  +  — 

V  2 

Y 

In  the  present  paper,  the  rate-dependent  isotropic  plasticity 
model  with  temperature-  or  humidity-dependent  material  prop¬ 
erties  is  introduced  into  the  FEM  and  solved  with  one  of  the 
commercial  FEM  software  packages,  such  as  ANSYS,  NAS- 
TRAN,  or  ABAQUS.  In  this  study  we  use  ANS  Y S  [1 8] ,  in  which 
the  CZM  and  the  contact  model  are  already  included. 


agglomerate  interface  after  a  long  period  of  operation  under  a 
cycled  humidity  or  temperature  load.  Interfacial  delamination 
can  be  modeled  by  traditional  fracture  mechanics  methods  such 
as  the  nodal  release  technique.  The  use  of  this  technique  requires 
one  to  select  the  crack  path  a  priori  so  that  double  nodes  (two 
nodes  sharing  the  same  spatial  position)  can  be  placed  on  it.  In 
order  to  avoid  this  constraint,  techniques  are  needed  that  directly 
introduce  one  kind  of  constitutive  law  for  describing  the  fracture 
mechanism  on  the  interfacial  surface,  which  is  taken  to  be  a 
phenomenological  mechanical  relation  between  the  traction  and 
displacement  jump  across  the  surface.  The  CZM  [19-25]  adopts 
softening  relationships  between  tractions  and  separations,  which 
in  turn  introduce  a  critical  fracture  energy  that  is  also  the  energy 
required  to  break  apart  the  interface  surfaces.  A  special  set  of 
interface  elements  is  used  to  represent  the  interface  surfaces  of 
the  materials,  and  the  CZM  is  used  to  characterize  the  interface 
surface  constitutive  behaviour  in  ANSYS  [18]. 

The  CZM  consists  of  a  constitutive  relation  between  the  trac¬ 
tion  acting  on  the  interface  and  the  corresponding  interfacial 
separation  8  (the  displacement  jump  across  the  interface,  which 
also  can  be  described  as  normal  separation  8n  and  shear  separa¬ 
tion  <St).  An  exponential  form  of  the  CZM,  originally  proposed 
by  Xu  and  Needleman  [20],  uses  a  surface  potential, 


cp(S)  =  eamax8n 


(5) 


where  0(8)  is  the  surface  potential,  e  the  Euler’s  number 
(2.7182...),  <rmax  the  maximum  normal  traction  at  the  inter¬ 
face,  8n  the  normal  separation  across  the  interface,  8t  the  shear 
separation  along  the  interface,  8n  the  normal  separation  where 
the  maximum  normal  traction  is  attained  with  8t  =  0,  and  8t  is  the 
shear  separation  where  the  maximum  shear  traction  is  attained 
at  8t  =  8t/ \fl.  So  the  traction  (on  both  the  normal  direction  and 
the  shear  direction)  can  be  defined  as 


Tn 


d<P 

dfn 


d& 

and  Tt  —  — . 

88t 


(6) 


From  Eqs.  (5)  and  (6),  the  normal  traction  of  the  interface  is 


2.3.  Material  models  for  interfaces  in  CLs 

There  are  interfaces  between  the  different  components  in 
CLs,  such  as  between  Nation  and  the  carbon/Pt  agglomerate, 
and  between  pore  and  solid  (as  indicated  in  Fig.  2).  The  rela¬ 
tionship  between  strain  and  stress  on  an  interface  is  obviously 
different  from  that  in  bulk  material.  So  in  our  simulation,  the 
CZM,  describing  the  behaviour  of  the  interface  between  Nation 
and  the  carbon/Pt  agglomerate,  and  the  frictional  contact  model, 
describing  the  behaviour  of  the  boundaries  between  pore  and 
solid,  are  proposed  to  solve  this  problem. 

2.3.1.  Cohesive  zone  model 

Fracture  or  delamination  along  an  interface  between  phases 
plays  a  major  role  in  limiting  the  durability  and  the  ductility 
of  multiphase  materials.  In  the  CLs  of  PEMFCs,  it  is  likely  that 
fracture  or  delamination  will  happen  along  the  Nafion-carbon/Pt 


and  the  shear  traction  is 


Fig.  3  shows  the  schematic  change  of  normal  and  shear  trac¬ 
tions  with  normal  and  shear  separations.  Both  the  normal  traction 
and  the  shear  traction  can  be  used  to  reveal  the  physical  graph 
of  interfacial  delamination  in  a  multiphase  material. 

In  our  simulation,  the  rate-dependent  properties  of  the  CZM 
are  not  explicit  when  compared  with  the  bulk  material.  However, 
when  the  parameters  in  the  CZM  are  connected  to  the  yield  stress 
under  different  loading  rates,  as  indicated  in  Section  2.2  and  Eq. 
(4),  the  CZM  can  also  be  considered  to  be  a  rate-dependent 
model. 
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Fig.  3.  Schematic  of  interfacial  traction  based  on  CZM:  (a)  normal  traction  and  (b)  shear  traction. 


It  should  be  mentioned  that  in  ANSYS  [18],  the  constitutive 
relation  for  each  cohesive  surface  is  taken  to  be  elastic  so  that 
any  energy  dissipation  associated  with  separation  is  ignored. 
Thus  the  failure  of  the  CZM  elements  is  manually  considered 
when  the  strain  in  the  CZM  exceeds  one  initial  strain  strength 
for  the  statistics  of  interface  failure  between  Nation  and  the  C/Pt 
agglomerate;  that  is,  for  one  CZM  element,  if 


max 


fSn  &t\ 

WStJ 


>  CO, 


(9) 


it  is  considered  to  be  a  failure  element.  8n  is  the  normal  separation 
across  the  interface,  8t  is  the  shear  separation  along  the  interface, 
and  co  is  the  coefficient  corresponding  to  the  maximum  strain 
of  interfaces.  Both  8n  and  8t  are  evaluated  from  the  average 
values  of  all  separations  at  Gauss  integration  points  in  the  CZM 
elements. 


nitude  across  their  interface  before  they  start  sliding.  This  state 
is  known  as  sticking.  The  Coulomb  friction  model  defines  an 
equivalent  shear  stress  r,  at  which  sliding  on  the  surface  begins 
as  a  fraction  of  the  contact  pressure  P, 

T  =  flP  +  P0, 

where  /x  is  the  coefficient  of  friction  (COF)  and  Po  specifies 
the  cohesion  sliding  resistance.  The  sticking/sliding  calculations 
determine  when  transitions  from  sticking  to  sliding  or  vice  versa 
occur. 

In  summary,  the  FEM,  with  the  support  of  the  CZM  and  the 
contact  model,  is  proposed  to  calculate  the  mechanical  response 
in  CLs.  The  underlying  mechanism  of  micro  structure  changes 
is  investigated  through  a  mechanical  analysis. 

3.  Parameters  in  material  models 


2.3.2.  Contact  model 

As  mentioned  in  Section  2.3.1,  the  CZM  is  suitable  only 
for  the  interface  between  different  solid  phases  in  CLs.  Contact 
models  must  be  used  to  describe  those  contact  phenomena  that 
emerge  with  environmental  changes.  Actually,  contact  models 
are  highly  nonlinear  and  require  significant  computer  resources 
to  solve  [26] .  Here,  a  contact  model  is  applied  on  the  boundary 
between  pore  and  solid  phases.  When  two  solid  domains,  con¬ 
sisting  of  perhaps  one  or  two  kinds  of  materials,  meet  because 
of  deformation,  the  contact  model  is  applicable.  On  the  newly 
formed  interface,  the  compress  force  can  be  transferred.  When 
there  is  tension  occurring  on  this  interface,  the  interface  sep¬ 
arates.  This  kind  of  contact  phenomenon  is  called  a  non-stick 
contact.  That  is,  for  one  contact  element,  if  on  >  0,  the  element 
is  considered  to  be  separated.  Here  on  is  the  normal  stress  along 
the  newly  formed  interface  between  two  domains. 

The  effect  of  friction  can  be  modelled  through  the  Coulomb 
friction  model  [26],  based  on  the  transition  between  the  stick¬ 
ing  and  sliding  states.  In  the  basic  Coulomb  friction  model,  two 
contacting  surfaces  can  carry  shear  stresses  up  to  a  certain  mag- 


Many  researchers  have  done  mechanical  experiments  on 
Nation  and  carbon  materials.  In  this  section,  some  parameters 
used  in  the  material  models  are  discussed  based  on  the  experi¬ 
mental  data  in  the  literature. 

3.1.  Mechanical  models  ofNafion  and  C/Pt  agglomerate 
materials 

One  of  the  most  representative  experiments  was  carried 
out  by  Tang  and  his  co-workers  [13].  In  their  experiment,  a 
MTS  Alliance  RT/5  material  testing  system  is  used  to  test  the 
mechanical  response  of  Nation  112  under  different  conditions 
by  adjusting  the  humidity  and  temperature.  One  of  their  exper¬ 
imental  results  (shown  in  Fig.  4)  gives  stress-strain  relations  at 
50%  RH  for  different  temperatures.  Similar  experiments  have 
been  done  by  other  researchers  on  Nation  115  [27-30],  Nation 
117  [29-35],  and  Nation  211  [30,36]. 

Based  on  the  experimental  data  in  the  literature,  a  piece- 
wise  linear  stress-versus-strain  curve  is  used  to  describe  the 
elastic-plastic  behaviour  ofNafion,  as  shown  in  Fig.  5.  From  the 
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Strain(mm  mm-1) 

Fig.  4.  Engineering  stress-strain  curves  of  tensile  tests  at  50%  RH  for  different 
temperatures  [13]. 

figure,  it  can  be  seen  that  six  parameters  are  needed  to  describe 
the  stress-strain  relationship  of  Nation :  the  elastic  Young’s  mod¬ 
ulus  E0,  the  plastic  modulus  E\,  the  static  yield  stress  apo,  the 
strength  as,  the  yield  strain  ep,  and  the  break  strain  es.  When 
the  Nation  stress  a  is  below  its  yield  stress,  Nation  behaves  as 
an  elastic  material;  when  it  is  between  its  yield  stress  and  its 


Fig.  5.  The  piecewise  linear  constitutive  relationship  of  Nation  adopted  in  the 
model. 

strength  stress,  Nation  behaves  as  a  plastic  material;  and  when 
it  exceeds  its  strength  stress,  Nation  exhibits  the  characteristics 
of  a  brittle  material.  Among  the  six  parameters,  only  four  are 
independent,  and  the  other  two  can  be  deducted.  For  simplic¬ 
ity,  here  we  chose  Eo,  tfpo,  crs>  and  ss  as  the  four  independent 
parameters.  The  Nation  yield  stress  and  strength  stress  can  then 


Fig.  6.  Material  parameters  as  functions  of  temperature  and  RH:  (a)  Young’s  modulus;  (b)  yield  stress;  (c)  strength;  (d)  break  strain.  The  points  are  experimental 
results  from  Refs.  [13,27,30,36]. 


F.  Rong  et  al.  /  Journal  of  Power  Sources  175  (2008)  699-711 


705 


Table  1 


Four  independent  parameters  in  constitutive  relationship 


a\ 

a2 

a3 

<34 

<35 

<36 

Error  (%) 

£o 

1.90E+00 

— 2.02E+01 

5.75E+04 

2.30E-02 

— 2.13E— 01 

3.87E+02 

0.435 

tfpO 

— 8.27E— 04 

— 9.60E— 03 

6.62E+01 

—  1.46E— 04 

—  1.29E— 03 

7.71E+00 

0.208 

1.33E-02 

— 6.41E— 02 

4.91E+02 

—  1.37E— 03 

8.00E-04 

1.82E+01 

0.570 

£s 

2.85E-05 

— 3.38E— 04 

3.47E+00 

7.77E-04 

— 2.42E— 03 

3.07E+00 

1.08 

be  expressed  as  follows: 


°p0  —  EoSp 

&s  —  VpO  —  E\(ss  —  £p) 


The  dependence  of  these  four  independent  parameters  on 
humidity  and  temperature  is  shown  in  Fig.  6.  Many  experiments 
have  been  done  to  study  the  effect  of  humidity  and  tempera¬ 
ture  on  these  parameters  [13,27-36].  Based  on  the  experimental 
data  ([13,27,30,36]),  the  dependencies  of  the  four  independent 
parameters  on  humidity  and  temperature  are  given  as 


Xr 


/RH 

VRH^ 


2 

+  Cl2 


+  a 3 


X 


(10) 


Fig.  7.  CME  value  as  a  function  of  T and  RH.  The  points  are  experimental  results 
from  Refs.  [13,27]. 


where  X  can  be  Eq,  crpo,  crs  and  £s;  at  (i=  1-6)  are  the  fitted 
parameters  from  the  experimental  data;  RH  and  T  are  the  current 
relative  humidity  and  temperature,  respectively;  and  RHr,  Tr, 
and  Xv  are  corresponding  reference  value  (here  the  reference 
values  are  input  as  the  values  at  30%  RH  and  25  °C).  The  fitting 
parameters  and  relative  errors  are  given  in  Table  1,  while  Fig.  6 
shows  the  results. 

In  addition,  the  Poisson  ratio  v  of  Nafion  can  be  estimated 
from  the  experimental  data.  Based  on  the  work  of  Tang  [13],  v  is 
not  sensitive  to  the  variation  of  temperature  and  humidity  and  is 
normally  taken  to  be  a  constant  (v  =  0.25)  [13,37].  The  swelling 
of  Nafion  that  occurs  under  different  humidity  levels  and  the 
expansion  of  Nafion  at  different  temperatures  can  be  character¬ 
ized  by  the  CME  and  CTE,  respectively.  Because  the  deviation 
of  the  CTE  is  relatively  small  compared  to  the  deviation  of  the 
CME  when  RH  and  temperature  change,  the  CTE  is  also  con¬ 
sidered  to  be  a  constant  in  this  paper  and  taken  as  1.23E— 4/°C 
from  Ref.  [37].  The  expression  of  the  CME  as  a  function  of  RH 
and  T  is  adapted  from  Ref.  [13]  and  [27],  based  on  Eq.  (10)  and 
shown  in  Fig.  7.  The  fitting  parameters  and  relative  errors  are 
given  in  Table  2,  in  which  the  reference  value  of  /3V  is  selected 
as  the  CME  at  30%  RH  and  25  °C. 

We  cannot  directly  obtain  the  rate-dependent  parameters  for 
Nafion  due  to  the  lack  of  experimental  data.  However,  many 
experiments  have  been  done  on  the  viscoplasticity  of  PTFE. 
Because  PTFE  is  the  backbone  of  Nafion,  it  is  reasonable  to 
assume  that  the  viscoplasticity  of  Nafion  is  similar  to  that  of 
PTFE.  Here  we  use  the  experimental  data  for  PTFE  to  get  the 
rate-dependent  parameters  of  Nafion  through  Eq.  (4).  The  exper¬ 
imental  results  for  PTFE,  taken  from  Ref.  [39],  are  given  in 


Fig.  8(a).  The  strain  rate  hardening  parameter  m  and  the  material 
viscosity  parameter  y  can  be  estimated  by  fitting  the  experimen¬ 
tal  results  with  Eq.  (4)  as 

m  =  0.230  ±  0.155;  y  =  0.0770  ±  0.0549, 

In  this  paper,  the  mean  values  are  chosen;  that  is,  hardening 
parameter  m  is  0.230  while  viscosity  parameter  y  is  0.0770. 
From  Fig.  8(b),  it  can  be  seen  that  the  ratio  of  the  rate-dependent 
yield  stress  to  the  static  yield  stress  increases  from  1.0  to  1.8 
when  the  loading  strain  rate  changes  from  0  to  1.0.  The  effect 
of  such  variation  on  the  microstructure  changes  of  CLs  will  be 
discussed  in  Section  4. 

Based  on  Ref.  [37],  the  C/Pt  agglomerate  is  stronger  than 
Nafion.  So  in  this  model,  the  C/Pt  agglomerate  is  considered  to 
be  an  elastic  material  and  there  is  no  failure  or  plasticity  accu¬ 
mulation  in  it.  The  Young’s  modulus  of  the  C/Pt  agglomerate  is 
taken  to  be  4.8  GPa  [37]. 


Table  2 

Parameters  of  humidity-  and  temperature-dependent  CME 


P 

a\ 

3.35E-10 

a2 

—  1.81E— 10 

<33 

1.17E-07 

<34 

4.87E-08 

<35 

—  1.14E— 07 

<36 

3.57E-04 

Error  (%) 

0.319 
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Fig.  8.  (a)  Strain-rate  dependence  of  PTFE  in  strain-controlled  tension  tests  using  different  strain  rates  [41];  (b)  the  ratio  between  rate-dependence  yield  stress  and 
static  yield  stress  of  Nation  changes  owing  to  strain  rate. 


3.2.  Interfacial  property  ofNafion  and  C/Pt  agglomerate 

The  interfacial  COF  pi  is  used  for  the  Coulomb  friction  model. 
The  COF  can  be  found  in  Ref.  [38]  and  is  shown  in  Fig.  9. 
Because  the  loading  rate  is  relatively  small,  only  friction  under 
low  loading  rate  is  considered.  So  pi  should  be  0.20  according 
to  Fig.  9. 

In  order  to  use  the  CZM  to  simulate  the  boundary  evolution 
between  Nation  and  the  C/Pt  agglomerate,  three  parameters  in 
Eq.  (5)  need  to  be  determined:  the  maximum  normal  traction  at 
the  interface  <7max,  the  normal  characteristic  separation  8n,  and 
the  shear  characteristic  separation  8t.  Empirical  relationships 
have  often  been  the  best  source  for  parameters  in  the  CZM,  and 
there  are  several  publications  on  determining  the  parameters  of 
the  interface  of  the  polymer  (such  as  epoxy  or  PTFE)  and  the 
carbon  nanotube  [39].  Based  on  these  works,  amax  is  close  to 
Eo/10,  where  Eq  is  the  elastic  modulus  of  the  material  (in  our 
case,  Nation).  The  normal  characteristic  separation  is  taken  to 
be  the  same  as  the  shear  characteristic  separation  [20] .  The  max¬ 
imum  shear  stress  is  close  to  crv/\/3  [18],  where  ap  is  the  yield 
stress  of  Nation.  The  shear  characteristic  separation  8t  is  calcu- 
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Fig.  9.  COF  of  Nafion  due  to  different  shearing  rates. 


lated  from  Eq.  (8).  Because  of  the  rate-dependent  characteristics 
of  yield  stress  (Eq.  (4)),  the  properties  of  the  CZM  can  be  also 
considered  rate-dependent. 

Eq.  (9)  is  used  to  determine  the  failure  of  the  CZM.  Because 
Nafion  is  weaker  than  the  C/Pt  agglomerate,  co  is  assumed  to  be 
the  failure  strain  of  Nafion,  which  can  be  calculated  from  Eqs. 
(4)  and  (10). 

It  is  should  be  mentioned  that  the  determination  of  some 
parameters  of  Nafion  and  the  C/Pt  agglomerate  is  very  difficult 
at  this  stage  and  we  have  borrowed  some  values  for  Nafion  from 
the  PTFE  material  in  order  to  carry  out  the  numerical  simulation 
and  to  analyze  qualitatively  the  mechanism  of  the  microstructure 
changes  in  the  CLs.  More  experimental  studies  or  molecular  sim¬ 
ulations  need  to  be  done  to  characterize  the  material  properties 
for  Nafion  and  the  C/Pt  mixture. 


4.  Numerical  analysis 

In  this  section,  we  simulate  the  microstructure  changes  in 
the  CLs  of  a  PEMFC  during  its  operation  by  applying  a  cycled 
humidity  and  temperature  load.  As  we  pointed  out  in  Section 
1,  the  microstructure  changes  might  result  from  chemical  prop¬ 
erty  variation  and/or  from  mechanical  property  changes.  As  a 
first  step,  in  this  paper  we  focus  only  on  the  changes  caused  by 
mechanical  property  variation. 

The  aim  of  the  simulation  is  to  understand  the  basic  mecha¬ 
nism  of  microstructure  changes  and  their  effects  on  performance 
based  on  the  model  proposed  in  Section  2.  To  verify  the 
CZM  and  the  contact  model,  one  microstructure  representa¬ 
tion  of  a  three-phase  microstructure,  including  Nafion,  C/Pt 
agglomerate  and  pore,  is  implemented.  The  relationship  between 
the  electrochemical  performance  of  fuel  cells  and  the  CL 
microstructure  is  a  complex  one,  with  factors  including  sur¬ 
face  area,  boundary  length  of  the  three  phases,  and  phase 
connectivity/tortuosity  playing  key  roles  [40] .  Thus  the  model 
simulates,  with  the  micro  structure  representative,  the  chang¬ 
ing  of  phase  connectivity  driven  by  the  cycled  humidity  and 
temperature. 
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Fig.  10.  (a)  A  typical  microstructure  representation  of  a  CL,  where  the  pink  domain  (domain  2)  represents  the  C/Pt  agglomerate  and  the  cyan  domain  (domain  1) 
represents  Nafion;  (b)  the  CZM  (red  elements)  and  contact  elements  (purple  line)  in  the  model.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the 
reader  is  referred  to  the  web  version  of  the  article.) 


4.1.  Computational  domains,  boundary  conditions,  and 
cycled  RH/T  loads 

A  typical  microstructure  representation  of  CLs  is  illustrated 
in  Fig.  10(a)  with  all  three  phases  present.  Microstructural  obser¬ 
vations  [5]  attest  that  the  contact  zones  between  the  solid  phases 
may  contribute  partially  or  completely  to  the  voids  in  the  whole 
CL.  In  Fig.  10(a),  Nafion  is  represented  by  the  cyan  domain 
(domain  1)  and  has  the  material  properties  described  in  Sec¬ 
tion  3.1.  The  C/Pt  agglomerate  (domain  2)  is  the  pink  domain 


and  has  elastic  material  properties.  Eight-node  structural  plane 
elements  (PLANE183)  are  used  to  make  up  the  domains  of  the 
bulk  material.  To  simulate  the  interface  of  Nafion  and  the  C/Pt 
agglomerate,  the  CZM  is  applied  to  the  boundary  elements, 
as  indicated  in  Fig.  10(b).  When  a  PEMFC  is  started  up  or 
shutdown,  the  humidity  and  temperature  in  the  CLs  changes 
periodically,  and  the  Nafion  swells  or  expands  accordingly. 

In  addition,  the  contact  model  is  applied,  as  shown  in 
Fig.  10(b).  Three-node  surface-to-surface  contact  elements 
(CONTA172)  are  applied  at  the  interfaces  between  solid  and 


Fig.  11.  Three  FEM  meshes  with  different  element  densities:  (a)  mesh  1  (865  elements);  (b)  mesh  2  (1739  elements);  (c)  mesh  3  (3824  elements). 
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pore,  while  the  boundary  between  different  solid  components  is 
meshed  by  six-node  cohesive  zone  elements  (INTER203).  In  our 
simulation,  the  plane  strain  condition  is  assumed.  The  symmet¬ 
ric  boundary  condition  is  applied  on  the  left,  right,  and  bottom 
boundaries  of  the  whole  computational  domain.  The  symmetric 
boundary  condition  is  also  used  on  the  top  boundary  of  Nation 
(domain  1),  while  the  free  boundary  condition  is  used  on  the  top 
of  the  C/Pt  agglomerate  (domain  2).  ^ 

The  temperature  and  humidity  profiles  come  from  Ref.  [12]. 

We  take  the  driving  force  for  the  microstructure  changes  in 
the  CLs  to  be  the  cycled  change  of  RH  and  T,  as  shown  in 
Fig.  1.  For  simplicity,  the  periods  of  humidity  and  temperature 
are  equal  (/rh  =  /r).  The  relative  humidity  starts  from  t-  0  while 
the  temperature  begins  after  a  2-s  delay  (to  =  2  s),  as  shown  in 
Fig.  1(b)  and  (c).  The  maximum  and  minimum  values  of  the 
RH  are  RHmax  =  90%  and  RHmin  =  30%,  respectively.  The  max¬ 
imum  and  minimum  values  of  temperature  are  Tmax  =  85  °C  and 
Tinin  =  25  °C,  respectively. 

4.2.  Mesh  dependency  of  the  model  and  phase  connectivity 
evolution 

The  CZMs  are  surface  elements  (line  elements  in  2D).  The 
kinematics  of  CZM  allows  for  displacement  jump  discontinu¬ 
ities  across  them.  And  the  elements  incorporate  a  constitutive 
law,  which  relates  the  displacement  discontinuity  to  the  cohesive 
force  to  be  transmitted  across  the  element.  Although  some¬ 
what  successful  in  yielding  numerical  results  that  correlate 
sufficiently  well  with  experimental  results,  this  method  still  con¬ 
strains  the  crack  path  to  evolve  only  along  the  interelement 
boundaries  of  the  underlying  FEM  grid.  Therefore,  the  present 
numerical  model  using  the  CZM  may  require  mesh  sensitivity 
studies  to  verify  that  the  underlying  FEM  grid  is  not  interfering 
with  the  physics  of  the  problem.  Fig.  1 1  shows  the  three  meshes, 
corresponding  to  different  numbers  of  elements  (865,  1739,  and 
3824,  respectively),  used  to  study  the  mesh  dependence  of  the 
present  model. 

In  addition,  although  a  comprehensive  analysis  of  the 
micro  structure  changing  performance  correlation  is  beyond  the 
scope  of  this  paper,  the  preliminary  analysis  below  demonstrates 
the  use  of  the  proposed  model  to  extract  key  microstructure  fea¬ 
tures.  Phase  connectivity  is  the  interface  between  Nation  and  the 
C/Pt  agglomerate  and  is  the  key  factor  forjudging  the  tortuosity 
of  protonic  and  electronic  transport  paths.  Thus  it  is  widely  con¬ 
sidered  in  determining  effective  contact  area,  effective  chemical 
reaction  rate,  and  percolation  limits.  From  the  proposed  model, 
the  connection  between  Nation  and  the  C/Pt  agglomerate  can 
be  determined  by  simply  calculating  the  contact  length  between 
them.  The  phase  connectivity,  that  is  the  contact  length  between 
different  solid  components  in  CFs,  can  be  measured  after  certain 
cycles  based  on  the  simulation,  as  shown  in  Fig.  12.  In  all  simu¬ 
lations,  the  periods  of  humidity  and  temperature  cycles  are  400  s. 

For  all  cases  in  the  present  simulations,  no  numerical  diver¬ 
gence  was  encountered.  Although  there  are  indeed  some  mesh 
sensitivities  in  the  present  method,  the  same  trend  can  be  seen 
in  all  three  simulations.  As  shown  in  Fig.  12,  for  the  simulation 
based  on  mesh  (a),  the  phase  connectivity  decreased  2.9%  of  the 


Fig.  12.  Decrease  of  connectivity  between  Nation  and  C/Pt  agglomerate 
(to  =tr  =  400  s)  based  on  three  different  mesh  densities. 

whole  initial  phase  connectivity  after  50  cycles.  For  mesh  (b),  the 
connectivity  decreased  3.2%,  while  the  connectivity  decreased 
3.6%  for  mesh  (c).  This  shared  phenomenon  cannot  be  caused 
by  numerical  error  accumulation  because,  based  on  the  conver¬ 
gence  of  the  FEM,  denser  FEM  grids  may  converge  to  the  true 
solution  rather  than  to  more  errors.  The  physics  of  the  prob¬ 
lem,  which  is  the  decrease  of  phase  connectivity  after  certain 
duty  cycles,  is  not  affected  by  the  mesh  sensitivity  of  the  model. 
Therefore,  the  convergence  and  stability  of  the  proposed  model 
can  be  proven  from  these  mesh  dependency  studies. 

Moreover,  this  kind  of  connectivity  between  Nation  and  the 
C/Pt  agglomerate  should  be  proportional  to  the  effective  electro¬ 
chemical  reaction  surface  in  the  CFs  based  on  the  three-phases 
surface  model.  In  addition,  this  decrease  indicates  the  possibil¬ 
ity  of  agglomerate  growth.  The  reduced  connectivity  induced 
by  the  cycled  humidity  and  temperature  change  could  indicate 
a  decrease  in  performance. 

4.3.  Basic  understanding  of  microstructure  changes  based 
on  parametric  study 

Because  there  are  so  few  mechanical  experiments  on  Nation 
and  the  interface  between  Nation  and  the  C/Pt  agglomerate 
and  the  proposed  model  is  hard  to  validate  by  experiment,  it 
is  necessary  to  study  the  effects  of  different  parameters  on  the 
microstructure  changes.  It  should  be  emphasized  that  our  pur¬ 
pose  is  to  analyze  qualitatively  the  mechanism  induced  by  the 
duty  cycles  based  on  the  parametric  study.  However,  this  para¬ 
metric  study  cannot  be  a  substitute  for  experimental  validation. 
On  the  contrary,  experimental  validation  is  urgently  needed  to 
confirm  our  understanding. 

Based  on  the  simulation  of  mesh  (b),  Fig.  13  shows  the  evo¬ 
lution  of  both  the  interfacial  delamination  fracture  energy  at  the 
boundary  between  different  components  and  the  plastic  Mises 
strain  in  Nation.  According  to  Fig.  13(a),  when  the  humidity  and 
temperature  begin  to  increase,  the  Nation  swells  and  the  delam¬ 
ination  fracture  energy  in  the  CZM  accumulates.  When  Eq.  (9) 
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Fig.  13.  (a)  The  delamination  fracture  energy  pattern  in  CZM  at  32  s;  (b)  the  delamination  fracture  energy  pattern  in  CZM  at  116  s;  (c)  the  delamination  fracture 
energy  pattern  in  CZM  at  2432  s  ((a)-(c)  share  the  same  legend);  (d)  the  plastic  Mises  strain  pattern  in  Nation  at  1 16  s.  There  is  no  plastic  strain  accumulated  in  Nation 
at  32  s;  (e)  the  plastic  Mises  strain  pattern  in  Nation  at  2432  s.  There  is  crack  initiation  in  Nation  close  to  the  top  of  the  boundary  between  the  different  components 
at  the  next  step  (7rh  =tr  =  400  s). 


is  satisfied  as  the  result  of  fracture  energy  accumulation  in  one 
CZM  element,  that  element  is  considered  to  be  the  failure  one. 
So  as  shown  in  Fig.  13(b),  the  element  at  the  triple  boundary 
of  Nafion,  C/Pt  agglomerate,  and  pore  failed  at  116  s.  Mean¬ 
while,  the  plastic  strain  accumulates  in  the  Nafion.  Because  the 
energy  in  the  Nafion  dissipates  when  the  CZM  element  fails  at 
116  s,  delamination  between  Nafion  and  the  C/Pt  agglomerate 
occurs.  At  the  same  time,  the  plastic  energy  is  also  accumulat¬ 
ing  in  the  Nafion.  There  is  no  plastic  strain  in  the  domain  of  the 
Nafion  at  32  s,  because  it  is  just  the  beginning  of  the  simula¬ 
tion.  Therefore  the  pattern  is  not  shown  in  Fig.  13.  The  plastic 
Mises  strain  pattern  at  116  s  is  shown  in  Fig.  13(d).  At  2432  s, 
the  strain  in  the  Nafion  at  the  top  of  the  interface  between  the 
different  solid  components  exceeds  its  strength  (as  indicated  in 
Fig.  13(e)).  Thus  there  is  crack  initiation  in  the  Nafion  from  that 
element. 

Therefore,  as  the  cycled  humidity  and  temperature  change, 
the  plastic  strain  accumulates  in  the  Nafion,  resulting  in  fatigue 
(Fig.  13(d)  and  (e)).  If  there  is  no  other  mechanism  for  the 
microstructure  changes,  this  kind  of  plasticity  accumulation 
results  in  crack  initiation  and  propagation  in  the  Nafion.  How¬ 
ever,  with  the  delamination  fracture  energy  accumulating  (in 
the  CZM  elements,  Fig.  13(a,  b,  c))  on  the  interface  between 
Nafion  and  the  C/Pt  agglomerate,  there  is  another  possible  reason 


for  microstructure  changes.  Thus  different  periods  of  humidity 
and  temperature  result  in  different  accumulations  of  residual 
stress  in  the  Nafion  as  well  as  in  different  energy  accumula¬ 
tions  for  delamination.  Crack  initiation  means  that  the  strain  has 
exceeded  the  break  strain  of  the  first  element  (whichever  Nafion 
element  or  CZM  element  it  is).  For  Nafion,  the  break  strain  fol¬ 
lows  Fig.  6(d)  and  the  failure  element  of  the  CZM  is  determined 
through  Eq.  (9).  So  the  competition  between  crack  initiation 
and  delamination  plays  an  important  role  in  the  micro  structure 
changes. 

Because  of  the  viscoplasticity  of  Nafion,  the  strain  rate  will 
affect  the  yield  behaviour  of  the  Nafion.  From  Fig.  8(b)  it  can 
be  seen  that  when  the  hardening  parameter  is  0.27,  the  yield 
stress  of  the  Nafion  increases  more  than  two  times  compared  to 
the  quasi-static  one.  So  different  cycle  times,  which  also  indicate 
the  strain  rate  of  materials,  influence  the  microstructure  changes. 

In  our  simulations,  four  different  periods  of  humidity  and 
temperature  are  used,  /rh  =  tj  =  100  s,  200  s,  300  s,  400  s.  Fig.  14 
shows  the  initiation  time  of  delamination  between  Nafion  and 
the  C/Pt  agglomerate  as  well  as  the  initiation  time  of  cracking 
in  the  Nafion  in  these  four  periods.  It  can  be  seen  that  different 
periods  result  in  different  initiation  times. 

The  shorter  the  cycle  period  of  temperature  and  humidity  is, 
the  earlier  delamination  occurs.  This  means  that  if  the  start-up 
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Fig.  14.  The  effect  of  different  periods  of  humidity  and  temperature  on  the  ini¬ 
tiation  time  of  delamination  between  Nation  and  C/Pt  agglomerate  (t&,  denoted 
by  the  square  dotted  line)  and  on  the  initiation  time  of  cracking  in  Nation  (tc i, 
denoted  by  the  circle  dotted  line). 

or  shutdown  of  PEMFCs  is  too  frequent,  delamination  between 
Nation  and  the  C/Pt  agglomerate  will  occur  more  easily.  On  the 
other  hand,  crack  initiation  in  Nation  has  a  completely  opposite 
tendency;  that  is,  a  shorter  cycle  period  leads  to  a  longer  crack 
initiation  time.  This  can  be  explained  as  follows.  When  the  cycle 
period  is  shorter,  the  strain  rate  of  the  Nation  is  greater  and  the 
corresponding  yield  stress  is  larger  too.  As  a  result,  it  is  relatively 
difficult  for  the  Nation  to  yield  and  thus  to  fail.  So  under  shorter 
periods  of  humidity  and  temperature,  Nation  is  more  resistant 
to  cracking  than  under  longer  periods.  Thus  there  is  competition 
between  Nation’s  yield-failure  and  Nafion-C/Pt  agglomerate 
delamination  because  of  the  viscoplasticity  of  Nation. 

Simulation  also  shows  that  cracks  always  start  from  delam¬ 
ination  between  Nation  and  the  C/Pt  agglomerate;  that  is,  from 
the  breakage  of  the  CZM  elements,  as  shown  in  Fig.  14.  This 
is  because,  based  on  the  material  properties  given  in  Section  3, 
the  CZM  is  relatively  weaker  than  the  Nation  material  under 
quasi- static  loading. 

In  summary,  there  is  competition  between  the  delamination 
energy  accumulation  (at  the  interfaces  between  different  solid 
components)  and  plastic  residual  strain  accumulation  (in  the 
domain  of  the  Nation).  This  competition  is  affected  by  the  period 
of  temperature  and  humidity  cycles  and  plays  an  important  role 
in  microstructure  changes,  resulting  in  either  crack  initiation 
in  the  Nation  or  delamination  between  Nation  and  the  C/Pt 
agglomerate. 

5.  Conclusions 

An  FEM,  with  the  support  of  the  CZM  and  fractional  con¬ 
tact  model,  is  proposed  to  simulate  the  microstructure  changes 
in  CFs  due  to  hydrothermal  duty  cycles  of  PEMFCs.  Simula¬ 
tions  are  undertaken  through  the  commercial  software  ANSYS 
by  using  a  user  subroutine.  Using  the  humidity-,  temperature-, 
and  rate-dependent  plastic  material  properties,  the  representa¬ 
tion  of  the  CF  microstructure  is  simulated.  Using  three  different 
mesh  densities,  it  is  found  that  although  this  method  has  little 


mesh  dependency,  a  shared  physical  phenomenon  can  be  seen  in 
the  results.  Microstructure  changes  do  occur  after  certain  duty 
cycles  as  the  result  of  competition  between  the  plasticity  energy 
accumulation  in  the  Nafion  and  the  delamination  energy  accu¬ 
mulation  at  the  interface  between  different  solid  components. 

Moreover,  based  on  the  parametrical  study,  we  found  that 
frequent  start-up  and  shutdown  of  fuel  cells  -  that  is,  shorter 
periods  in  RH  and  T  cycles  -  lead  to  an  earlier  initiation  time  of 
delamination  between  the  Nafion  and  C/Pt  agglomerate,  but  a 
later  initiation  time  of  fracture  in  the  Nafion.  Finally,  the  length 
of  the  interface  between  Nafion  and  the  C/Pt  agglomerate  is 
also  investigated  after  certain  hydrothermal  cycles.  It  is  shown 
that  the  length  of  the  connection  decreased,  which  may  indicate 
performance  degradation. 

From  the  simulation  results,  it  is  found  that  the  model 
itself  recurs  what  the  experiment  show  in  its  representation 
of  micro  structure  in  CFs.  However,  CF  microstructure  is 
naturally  random  and  heterogeneous.  Therefore,  more  practical 
simulations  of  CF  microstructure  will  help  us  to  understand  the 
real  problem  more  accurately.  Work  is  now  ongoing  to  construct 
the  microstructure  based  on  information  from  practical  CFs  and 
to  analyze  the  microstructure  changes  due  to  the  duty  cycles  of 
fuel  cells. 
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